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ABSTRACT 

The  effective  design  of  body  armor  requires  the  use 
of  biomimicking  materials  that  accurately  and  robustly 
reproduce  the  characteristic  mechanical  behavior  of 
human  tissues.  The  mechanical  performance  of  the  soft 
tissues  is  determined  by  their  viscoelastic  properties  and 
for  this  reason  the  optimal  soft  tissue  surrogate  must 
viscoelastically  mimic  the  original  tissue.  One  traditional 
tissue  surrogate  material  is  a  water-based  gelatin  gel  that 
has  poor  environmental  stability  and  limited  tailorability. 
Gel  systems  based  on  self-assembled,  amphiphilic  ABA 
triblock  copolymers  have  similar  properties  to  the  gelatin 
system,  forming  the  stable,  spatially  extended  networks 
with  a  tunable  viscoelastic  behavior.  In  the  present  work, 
we  have  developed  a  systematic  way  to  evaluate  the 
viscoelastic  properties  of  the  biomimetic  gels  to  control 
their  mechanical  properties.  We  have  validated  our 
approach  using  experimental  data  and  demonstrated  that 
our  simulation  results  are  in  good  qualitative  agreement 
with  the  experimental  data.  The  structure-property 
relationships  found  from  the  mesoscale  simulation  allows 
us  to  control  mechanical  properties  of  these  gels  in  order 
to  mimic  a  wide  range  of  soft  tissues.  The  viscoelastic 
properties  have  been  calculated  employing  a  non¬ 
equilibrium  oscillatory  shear  technique  used  with  the 
Dissipative  Particle  Dynamics  method  (DPD). 

1.  INTRODUCTION 

The  mechanical  response  of  living  tissue  is  crucial  to 
identifying  the  potential  injuries  related  to  impact  events. 
The  U.S.  Army  still  routinely  uses  gelatin  as  a 
standardized  medium  to  test  the  terminal  performance  of 
firearms  ammunition.  Ballistic  gelatin  is  a 
thermoreversible  material  composed  of  denatured 
collagen  and  when  cooled  lower  than  30°  C,  it  forms  a  gel 
consisting  of  the  rigid  triple  helix  zones  interbridged  by 
flexible  coils.  The  ballistic  gelatin  approximately 
simulates  the  density  and  viscosity  of  human  soft  tissues 
and  provides  similar  performance  for  most  ballistics 
testing.  The  elastic  modulus  of  ballistic  gelatin  is 
approximately  equal  to  100-150  kPa  at  a  temperature  near 
10  °C,  and  is  comparable  with  the  elastic  modulus  of  soft 
tissue,  which  varies  from  25  to  300  kPa.  Low  cost  and 


availability  are  also  major  advantages  of  the  ballistic 
gelatin. 

However,  the  use  of  ballistic  gelatin  has  several 
significant  disadvantages.  Variable  molecular-weight 
distribution  of  naturally  produced  gelatin  and  procedural 
differences  can  lead  to  inconsistent  mechanical  properties 
and  impede  data  analysis.  The  evaporation  of  water  can 
also  alter  mechanical  properties  of  the  gelatin.  Finally,  as 
mentioned,  the  properties  of  ballistic  gelatin  only 
approximate  soft  tissue  over  a  narrow  temperature  range 
around  10  °C,  rapidly  losing  its  solid-like  consistency  as 
temperature  increases  toward  its  gel  temperature,  30  °C, 
at  which  point  it  becomes  a  liquid  (Juliano  et  al.,  2006). 

Taking  into  account  all  these  factors,  the  application 
of  novel  polymer  surrogate  systems  for  impact  studies  is 
needed  to  improve  the  repeatability  and  reproducibility  of 
testing  while  providing  robust  model  validation 
capabilities.  Polymer  gels  based  on  blends  of  synthetic 
copolymers  in  a  solvent  exhibit  improved  thermal  and 
environmental  stability  and  possess  tunable  viscoelastic 
properties  over  a  broad  range  of  strain  rates. 

ABA  triblock  copolymers  are  one  example  of  a 
polymer  gel  which  has  been  widely  used  in  many 
practical  applications.  ABA  triblock  copolymers  undergo 
microphase  separation  in  midblock-selective  solvent  and 
form  a  spatially  extended  network.  The  microstructure  of 
this  network  depends  on  several  factors  such  as 
temperature,  copolymer  concentration,  block  architecture, 
and  relative  length  of  the  hydrophobic  and  hydrophilic 
blocks.  The  polymer  network  consists  of  physical  cross¬ 
links  comprising  of  the  solvent-incompatible  endblocks 
(micelles)  connected  by  the  midblocks  (bridges).  The 
micelles  are  surrounded  by  a  corona  consisting  of 
midblocks  (loops)  formed  as  a  result  of  close  association 
of  the  endblocks  (Fig.  1).  The  bridges  act  as  active  chains 
of  the  network  and  are  treated  as  elastic  springs,  where 
loops  do  not  contribute  to  the  gel  network  elasticity 
(Drozdov  et  al.,  2004;  Tanaka,  2002).  For  that  reason, 
bridge-to-loop  ratio  is  one  of  the  important  factors 
defining  the  mechanical  properties  of  a  gel  and  this  ratio 
depends  on  several  factors  such  as  copolymer 
concentration  and  the  relative  length  of  the  blocks.  The 
bridge  to  loop  ratio  increases  with  relative  length  of  the 
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midblock  with  respect  to  the  length  of  the  endblocks  and  cNt 

copolymer  concentration  (Shen  et  al.,  2007;  Szczubialka  ^  ~  M 

et  al.,  2000). 
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In  this  study,  we  attempt  to  create,  through 
computational  modeling,  a  predictive  tool  that  answers 
questions  regarding  the  impact  of  chemical  and 
compositional  variables  on  the  viscoelastic  properties. 
Understanding  such  structure-property  relationships  will 
enable  the  development  of  biomimetic  gels  that  mimic  a 
wide  range  of  soft  tissues  over  a  wide  range  of  strain  rate. 


Fig.  1.  Schematic  representation  of  a  polymer 
network  formed  from  ABA  triblock  copolymer. 


where  M  =  x  +  2  and  C  is  a  volume  fraction  of  the 
triblock  in  the  solution.  It  is  important  to  mention  that  we 
model  an  endblock  of  any  size  as  a  single  bead  to  avoid 
its  collapse  to  a  smaller  size  before  the  association  with 
the  other  endblocks  takes  place. 

The  physical  size  of  the  interaction  radius  for  the 
DPD  particle,  Rc  depends  on  the  molar  volume  of  the 
endblock  and  may  be  estimated  from 


(2) 


where  p  is  the  DPD  bead  density,  vA  is  the  molar  volume 
of  the  endblock,  and  NA  is  Avogadro’s  Number,  x  is 
evaluated  from  the  volume  fracture  of  the  end  block  in  the 
dry  copolymer,  fA 

x=2(1-^).  (3) 
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2.  MODEL  AND  SIMULATION  METHODS 

Generally  speaking,  atomistic  molecular  dynamics  is 
suitable  up  to  100  nanometers  and  several  nanoseconds, 
which  makes  it  difficult  to  perform  simulations  of  large 
polymeric  system  on  comparable  timescales.  In  this  study, 
we  have  exploited  the  Dissipative  Particle  Dynamics 
(DPD)  simulation  technique  which  consists  of  coarse¬ 
grained  particles  that  represent  molecular  clusters  rather 
than  individual  atoms.  The  DPD  method  is  a  mesoscale 
simulation  technique  that  operates  at  time  and  length 
scales  larger  than  those  of  traditional  molecular  dynamics, 
but  for  situations  that  are  still  inaccessible  at  the 
continuum  level. 

2.1  Coarse-Grained  Model  of  AiBxAxTriblock 
Copolymer 

The  three-dimensional  system  consists  of  Nt  DPD 
particles  of  triblock  copolymers  plus  solvent  confined  in  a 
cubic  domain  with  periodic  boundaries  in  all  directions. 
The  system  has  N  chains  of  AiBxAi  triblock  copolymer 
composed  of  particles  of  types  “A”  and  “B”  of  total 
length  M  immersed  in  a  midblock-selective  solvent  whose 
particles  are  of  type  “S.”  Each  DPD  bead  represents 
roughly  the  same  volume.  Each  chain  has  two  endblocks 
represented  by  a  single  bead  of  type  “A”  and  a  midblock 
composed  of  x  particles  of  type  “B”  (Fig.  2).  The  number 
of  triblock  chains  is  chosen  from 


Note  that,  because  the  endblock  of  any  length  is 
represented  by  one  DPD  bead,  x  corresponds  to  the 
relative  length  of  the  middlebock  in  the  triblock 
copolymer. 

The  poly(styrene-block-isoprene-block-styrene) 
copolymer  or  SIS  in  the  I-selective  solvent  has  been 
chosen  as  a  model  triblock  copolymer  for  this  study 
because  of  its  promising  properties  as  a  biomimetic  gel 
and  sufficient  amount  of  existing  experimental  structural 
and  rheological  data  to  validate  our  model. 


Fig.  2.  Schematic  representation  of  A1B3A1  triblock 
copolymer  mapped  on  DPD  model.  Poly(styrene-block- 
isoprene-block-styrene)  copolymer  is  chosen  as  an 
example. 
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2.2  Dissipative  Particle  Dynamics 

A  DPD  system  is  composed  of  soft  particles,  each 
representing  a  region  of  fluid,  which  moves  according  to 
Newton’s  equations  of  motion  continuously  in  space  and 
discreetly  in  time.  In  a  DPD  simulation,  the  polymer 
chain  is  modeled  as  a  collection  of  point  particles  that 
represent  lumps  of  the  chain  containing  several  segments. 
Force  acting  on  DPD  particles,  F  is  written  as  the  sum  of 
a  conservative  force  Fc,  dissipative  force  FD ,  and  random 
force  Fr : 

F  =  FC  +Fd  +Fr  (4) 

ij  y  ij  ij  v  7 


various  effects  (e.g.,  polymer  chain  length)  are 
manifested.  The  viscoelastic  properties  are  calculated 
using  a  non-equilibrium  oscillatory  shear  technique 
(Pryamitsyn  and  Ganesan,  2006),  which  entails  a 
simulation  with  an  additional  force  in  the  shear  direction 
along  with  time-dependent  Lees-Edwards  boundary 
conditions  (Lees  and  Edwards,  1972).  For  oscillatory 
shear  imposed  in  the  x-y  plane,  the  equation  of  motion  for 
the  particle  velocities  becomes 


dvix 


d2r(t) 

dt 2 


(7) 


The  oscillatory  strain  can  be  taken  as 
y(t)  =  Aa\[-cos(cQt)\,  where  Aa  and  co  are  the  chosen 
values  of  the  amplitude  and  frequency,  respectively. 


Fc  includes  a  soft  repulsion  force  FCr  acting  between 
two  particles  and  a  harmonic  spring  force  FCs  acting 
between  adjacent  particles  in  a  polymer  chain.  FCr  and  FCs 
are  given  by 


Ff  = 


and 


i — — 


AC  J  Aij 


r. 


for  r.  <  r 


for  r.  >  r 


(5) 

(6) 


respectively. 

In  Eq.  5  and  6,  ay  is  the  maximum  repulsion  between 
particle  i  and  particle  y,  ri-=ri—Y.,  Fj  “  |hj  |  ?  rc  r° 

cutoff  radius,  K  is  the  spring  constant,  and  r0  is  the 
equilibrium  spring  length.  The  evolution  of  DPD 
particles  in  time  t  is  governed  by  Newton's  equations  of 
motion.  For  a  more  detailed  description,  see  the  original 
papers  (Hoogerbrugge  and  Koelman,  1992;  Koelman  and 
Hoogerbrugge,  1993). 


2.4  Simulation  Parameters 

DPD  reduced  units  system  are  adopted  for  the 
convenient  expression  of  parameters  and  are  taken  as: 
length  in  Rc ,  energy  in  kBT/Rc ,  mass  in  particle  mass  m , 
and  time  in  t  =  Rc-sjm/kBT  (Groot  and  Rabone,  2001). 

Therefore,  the  physical  value  of  moduli  are  calculated 
from 


G  =  GDPDkBT/{Rc)\ 


(8) 


where  T  is  the  temperature  and  kB  is  Boltzmann's 
constant.  We  choose  m=kBT=Rc=  1,  time  step  At=0.02, 
and  overall  particle  density  p=3.  K  and  r0  in  Eq.  6  for  the 
triblock  are  1  and  3,  respectively. 

The  repulsive  parameters  for  like  particles  an  are 
determined  from  the  value  of  the  dimensionless 
compressibility  (Groot  and  Warren,  1997): 
k  ~  1  +  0.2 au  /  pT  ^  w]iere  £-i  js  gjven  by 


r1 


1 

nkBTkT 


(9) 


2.3  Viscoelastic  Property  Calculations 

When  a  linear  viscoelastic  material  is  subjected  to 
oscillatory  strains,  y  of  frequency  co ,  the  stress  response, 
a,  is  necessarily  cyclic  and  can  be  written,  as:  cj(co)  =  A  a 
(G'(co)sin(cot)+  G”(co)cos(cot)).  The  storage  modulus,  G\ 
indicates  the  material’s  ability  to  store  energy  and  the  loss 
modulus,  G”(co)  characterizes  the  amount  of  energy  lost 
through  the  viscous  process. 


and  where  kT  is  the  isothermal  compressibility  and  n  is  the 
number  density.  au  is  taken  as  25  for  water,  while  for 
organic  solvent  composed  of  larger-size  molecules  we 
estimated  a  a  to  be  50. 

The  repulsive  parameters  for  unlike  partciles  ay  are 
determined  according  to  the  linear  relation  with  Flory- 
Huggins  parameters  (Groot  and  Warren,  1997): 


Simulations  can  be  used  to  determine  which 
parameters  govern  the  viscoelastic  behavior  of  polymer 
systems,  delineating  the  regimes  and  frequencies  at  which 
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a,j  ~  au  +  3.27^. 


(10) 


where  %  parameters  are  mapped  from  the  real  polymers 
to  the  short  simulated  chains  preserving  %M  =  const . 

For  our  model  block  copolymer  system  SIS  in  an 
organic  solvent,  the  interaction  parameter  for  the 
polystyrene-polyisoprene,  Xps-pi  has  been  taken  from  the 
literature,  equal  to  0.086  for  30°  C  for  M  =  432. 
Interaction  parameter  between  endblock  and  solvent  is 
chosen  such  that  jM  >  50  to  ensure  good  segregation.  The 
maximum  repulsive  force  ay  is  calculated  according  to 
Eq.  10  and  given  in  Table  1. 


Table  1.  Repulsive  parameters  between  a  pair  of  DPD 
particles.  _ _ _ _ _ _ 


DPD  pair 

s-s 

A-A 

B-B 

A-B 

A-S 

B-S 

50.0 

50.0 

50.0 

66.0 

105.0 

50.0 

2.5  Methodology 

In  this  study,  we  have  investigated  the  role  of  relative 
block  size  and  concentration  of  triblock  copolymer  on  the 
structure  and  structure-related  mechanical  properties  of 
formed  polymeric  gels.  In  order  to  validate  our  method 
we  choose  in  B  selective  solvent  as  a  prototype  of 

SIS  of  molecular  weight  (7.2-35.8-7.2  kg/mol)  in  the  I- 
selective  solvent,  ?z-tetradecane  studied  by  Watanabe  et  al. 
(Watanabe  et  al.,  2000)  using  methodology  described  in 
Section  2.1.  We  choose  concentration  of  AiBgAi  0.2,  0.3, 
0.4,  and  0.5  wt,  which  correspond  to  the  volume  fraction, 
c  equal  to  0.17,  0.28,  0.39,  and  0.50.  The  triblock  would 
adopt  micellar  morphology  for  the  room  temperature  for 
this  concentration  range.  In  order  to  probe  effect  of 
relative  block  length,  we  consider  x  equal  to  3,  6,  9,  12, 
15,  and  18. 

We  have  performed  the  equilibrium  simulation  of  the 
system  consisted  of  randomly  generated  polymer  chains 
and  solvent  particles  slowly  adjusting  the  periodic  box 
size  to  the  desired  pressure  followed  by  a  constant 
temperature  annealing  simulation.  This  procedure  has 
been  repeated  six  times  verifying  that  the  system  pressure 
does  not  change  as  an  indication  of  an  equilibrium  state. 
Ten  parallel  simulation  runs  were  issued  for  each  set  of 
initial  conditions.  We  examined  the  microstructure  of 
each  equilibrium  state  to  determine  the  conformation  of 
each  midblock  in  the  gel  using  the  following  micelle- 
identifying  algorithm.  The  midblock  is  considered  to  be  in 
a  loop  conformation  if  both  endblocks  belong  to  the  same 
micelle.  If  the  midblock  connects  two  different  micelles  it 
is  considered  to  be  in  a  bridge  conformation.  An 
unassociated  endblock  can  also  form  a  dangling  tail 
conformation  for  the  midblock.  Visual  inspection  of  the 
microstructure  was  also  performed. 

Following  equilibration  of  the  structure,  a  series  of 
stress  trajectories  are  generated  by  imposing  the 


oscillatory  shear  conditions  with  0.1%  strain.  Since  the 
DPD  simulation  is  generally  able  to  access  only 
frequencies  in  the  megahertz  regime,  dynamic  moduli 
have  been  linear  extrapolated.  Note  that  numerical  errors 
tend  to  be  higher  at  low  frequencies  requiring  averaging 
over  more  stress  trajectories  to  minimize  these  errors, 
which  is  ineffective  and  leads  to  very  long  simulation 
times.  This  suggests  that  a  new  computational  approach 
to  estimating  the  mechanical  properties  of  ABA  triblock 
copolymer  is  needed. 

3.  RESULTS  AND  DISCUSSION 

3.1.  Validation  of  our  method. 

For  the  validation  purposes,  we  have  compared  the 
microstructure  and  the  mechanical  properties  of  our  gel  to 
experiment  (Watanabe  et  al.,  2000).  All  simulation  runs 
resulted  in  microphase  separation  which  led  to  a  three- 
dimensional  network  structure  of  a  micellar  gel  with  the 
middle  block  exhibiting  bridge  or  loop  conformations. 
Fig.  3.  shows  an  example  of  an  equilibrated  system. 


Fig.  3.  Equilibrated  structure  of  DPD  simulation  of 
A^Ai  copolymer.  Micelles  comprised  of  type  "A" 
particles  are  shown  in  blue.  The  bridge  and  loop 
conformations  of  the  middle  blocks  are  shown  as  red  and 
brown,  respectively.  Solvent  particles  are  not  shown. 

Fig.  4  compares  the  estimated  loop  fractions  of  the 
middle  block  from  the  DPD  simulation  (b),  with  those 
obtained  by  Watanabe  et  al.  by  measuring  the  dielectric 
losses  of  the  SIS  in  ^7-tetradecane  (a).  Simulation  results 
are  in  good  qualitative  agreement  with  the  experimental 
data. 
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0/  wt% 

Fig.  4.  Loop  fraction  as  a  function  of  copolymer 
concentration  (a)  for  SIIS  (7.2-17.9-17.9-7.2)/C14 
(Watanabe  et  al.,  2000)  (b)  from  mesoscale  simulation  of 
AiBgAi  triblock  copolymer. 

The  storage  modulus,  G '  obtained  by  our  simulations 
exhibits  quite  extensive  behavior  with  a  numerical  value 
close  to  a  plateau  for  accessible  DPD  frequencies.  We  do 
not  consider  the  loss  modulus  G"(co)  is  this  study  because 
its  value  is  prone  to  numerical  error  of  calculation.  The 
equilibrium  value  of  storage  modulus,  G  found  from  DPD 
simulation  has  been  mapped  back  to  its  physical  value 
using  Eq.  8.  Table  2  shows  the  simulated  and 
experimental  rheological  data  and  demonstrates  relatively 
accurate  representation  of  experimental  data  with  DPD 
simulation.  Note  that  because  the  I  blocks  in  the  SIS 
solutions  are  not  entangled  due  their  low  molecular 
weight  (Mw  ~  50000),  DPD  simulation  results  are  in  a 
good  agreement  with  experimental  data. 


Table  2.  Comparison  between  equilibrium  values  of 
shear  modulus  estimated  from  simulation,  Gsim  and 
obtained  from  experimental  data,  Gexp _ _ 


conc.,wt 

0.2 

0.3 

0.4 

0.5 

log(Gsim) 

3.29 

3.79 

3.95 

4.11 

log(GeXp) 

2.70 

3.40 

3.85 

4.20 

3.2  Governing  Equation  for  Elastic  Shear 
Modulus 

Gels  are  characterized  by  a  storage  modulus,  G '  that 
often  is  only  weakly  dependent  on  frequency,  co  over  a 
broad  range  and  is  several  order  of  magnitude  larger  than 
the  loss  modulus,  G"(co).  Thus,  elastic  response  appears 
to  be  a  reasonable  assumption,  such  that  the  equilibrium 
storage  modulus,  G  is  given  by 

G  =  PivhGh  +  viGi  +  V,G, )  (ll) 

where  p  is  copolymer  density,  vb,  V/  and  v*  are  number 
density  of  bridges,  tail  and  dangling  tails  and  Gb ,  G/  and 
Gt  are  the  elastic  modulus  of  a  single  bridge,  loop  or 
dangling  tail  respectively. 

For  the  strong  separation  condition  we  can  assume 
that  the  contribution  of  tails  is  negligible  and  that  all 
midblocks  have  either  bridge  or  loop  conformation.  Using 
Eq.  1  and  evaluating  number  of  bridges,  Nb  as 


Nb  =  N(j)b{x,c)  and  number  of  loops,  Ni  as 
we  get  the  equilibrium  elastic 

modulus  as: 

G  =  P  [G A  {x ,  c )  +  G,  (1  - 1  (x ,  c ))]  (12) 

x  +  2 

for  %M  >  50  ,  where  (f)b  (x,  c )  is  bridge  fraction,  that  is  a 

function  of  the  relative  block  size  and  copolymer 
concentration. 

3.3  Fraction  of  Bridges 

To  compare  effect  of  concentration  and  relative  block 
size  on  fraction  of  bridges,  (f)b  we  have  performed  a  set  of 

simulations  with  variable  concentration,  c  and  length  of 
the  midblock,  x.  We  present  our  results  in  Fig.  5.  From 

this  figure,  one  can  see  that  (j)b  increases  with  increasing 

of  the  copolymer  concentration.  (f)b  increases  with 

increasing  of  the  relative  length  of  the  middle  block  for 
the  shorter  middle  blocks  (x  <  9).  For  the  relatively  long 

middle  block  compared  with  the  endblock  (x  >  9),  (j)b  only 

weakly  depends  on  the  middle  block  length.  This 
behavior  is  especially  pronounced  for  low  concentration 

(< c  =  0.2),  where  (f)b  ~  0.4  for  all  x  >  9.  This  is  a  result  of 

the  lower  probability  for  the  endblocks  of  the  different 
copolymer  chains  in  a  dilute  triblock  solution  to  form  a 
network  since  close  association  is  minimized. 

The  compact  size  of  the  middle  block  compared  to 
the  distance  between  aggregates  leads  to  the  tendency  to 
form  loops.  For  example,  for  x  =  3,  all  middle  blocks  form 
loops  for  the  entire  concentration  range.  Increasing  the 
copolymer  concentration  increases  the  probability  of  each 

chain  developing  a  bridge  such  that  (/)b  reaches  a 

maximum  value,  ^max  for  the  the  pure  copolymer.  For  the 

reasonable  degrees  of  incompatibility 

XM  >  100  corresponding  to  the  relatively  long  midblock 

(x  >  9)  ^max  is  approximately  0.66.  This  value  is  in  good 
agreement  with  theoretical  work  (Jones  et  al.,  1996) 
having  ^max  equal  to  0.63. 

For  relatively  high  degree  of  incompatibility  (jM  > 
50),  ^is  only  a  function  of  concentration  and  triblock 

architecture.  However,  in  general,  (f)b  also  depends  on 

degree  of  incompatibility  or (j)b  =  This 

dependence  becomes  significant  for  endblocks  with  a  low 
degree  of  polymerization  in  a  non-highly  selective 
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solvent.  For  this  case,  the  number  of  unassociated 
endblocks  would  be  higher  and  both  (f)b  and  ^  would 

have  smaller  values.  To  address  this  situation  we  can 
change  the  repulsive  parameters  in  our  model  according 

to  Section  2.4  and  then  evaluate^  .  It  is  important  to 

mention  that  (pb  does  not  depend  on  the  absolute  length  of 
the  blocks  but  depends  on  the  relative  block  sizes. 

1 

0.8 

0.6 
_Q 

•e 

0.4 
0.2 
0 

0  5  10  15  20 

x 

Fig.  5.  Fraction  of  bridges  for  copolymer  as  a 

function  of  the  relative  block  size  at  the  various 
concentrations. 

3.4  Elastic  Shear  Modulus 

To  assess  the  effect  of  concentration  and  relative 
block  size  on  the  equilibrium  elastic  modulus,  G  we  have 
performed  oscillatory  shear  simulations  for  the  same  set 
of  concentration  and  length  of  the  middle  block  we  used 
for  calculation  of  number  of  bridges.  Results  of  shear 
simulations  for  AiBi5Ai  are  shown  in  Fig.  6. 
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Fig.  6.  Storage  modulus  for  A1B15A1  triblock 
copolymer  as  a  function  of  frequency  at  various 
concentrations. 


from  Eq.  12  setting  (f)b  equal  to  0.  Assuming  that  G/  does 

not  change  with  concentration  and  the  polymer  length  we 
found  that  G/  ~  0.1 14  in  DPD  units.  Having  this  value,  we 
are  able  to  find  from  Eq.  12  that  the  elastic  modulus  of  a 
single  bridge,  Gb  for  the  relative  long  middle  block  (x  >  9) 
does  not  depend  on  concentration  and  block  length  and 
equal  Gb  ~  1.717.  For  the  shorter  middle  block  (x  =  6)  Gb 
~  0.995  in  DPD  units.  From  these  values,  we  conclude 
that  contribution  of  the  loops  compared  to  the  bridges  to 
the  elastic  modulus  is  negligible. 

From  entropic  elasticity  theory  (Ferry,  1980),  the 
contribution  of  one  bridge  to  the  equilibrium  elastic 
modulus  is  estimated  as 


Gb  =  ^kBT , 


(13) 


where  ju  is  a  coefficient  of  order  of  unity.  Thus,  we  have 
found  that  ju  ~  1 .7 17  for  x  >  9. 

Mapping  our  results  to  physical  values  of  the  elastic 
modulus  for  the  ABA  triblock  copolymer  and  having 
specified  p= 3,  we  can  finally  write 


G  = 


1.717  RT 

va(x  +  2) 


c<pb{x,c,%M), 


(14) 


where  R  is  the  universal  gas  and  x  >  9.  In  Eq.  14  the 
volume  fraction,  c,  the  molar  volume  of  A  block,  vA  and 
degree  of  incompatibility  jM  characterize  the  physical 
triblock  copolymer  and  solvent.  The  relative  size  of  the 
midblock  with  respect  to  the  endblock,  x,  is  evaluated 
from  Eq.  3.  Corresponding  to  these  parameters,  the 

fraction  of  bridges  (j)b  (x,  c,  )  is  determined  from 

equilibrium  DPD  simulation.  Eq.  14  (solid  lines  in  Fig.  7) 
is  applicable  for  the  relative  long  midblock  compared  to 
the  endblock  size  (x  >  9).  For  x  <  9,  we  found  the 
coefficient,  p,  being  less  1.717,  for  example 
approximately  1.0  for  x  =  6. 

Because  (f)b  increases  with  increasing  relative  length 

of  the  middle  block  for  x  <  9  and  only  weakly  depends  on 
the  middle  block  length  for  x  >  9,  the  elastic  modulus  for 
each  concentration  has  a  maximum  value  and  then 
decreases  with  increasing  x  (Fig.  7).  According  to  Eq.  14, 
the  elastic  modulus  increases  with  concentration  directly 

and  indirectly  since  (f)b  is  also  increasing  with 
concentration. 


Following  the  procedure  described  earlier,  we 
obtained  G  for  various  concentrations  and  lengths  of  the 
midblocks,  and  results  are  presented  in  Fig.  7.  Since  the 
midblock  of  A^Ai  forms  only  loops  we  can  estimate  the 
elastic  modulus  of  a  single  loop  to  the  elastic  modulus,  G/ 


The  main  advantage  of  our  model  is  that  it  can 
estimate  the  elastic  response  of  ABA  triblock  copolymers 
of  low  molecular  weight  in  a  midblock-selective  solvent 
directly  from  equilibrium  DPD  simulations 
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evaluating  (j)b  (x,  c,  ) .  Since  non-equilibrium  shear 

simulation  is  not  necessary,  this  approach  greatly  reduces 
the  simulation  time. 


Fig.  7.  Elastic  modulus  for  AiBxAi  as  a  function  of 
the  relative  block  size  at  various  concentrations. 
Simulation  data  are  represented  with  points  and  solid  lines 
are  fit  to  Eq.  14. Vertical  lines  show  the  applicable  part  of 
our  model. 


For  the  second  example,  we  have  considered  a  blend 
of  diblock  and  triblock  copolymers  of  different  relative 
block  size  (one  gel  contain  15  %  by  mass  polysterene,  PS- 
15  and  the  other  gel  contains  30  %  by  mass  polysterene, 
PS-30)  studied  by  Juliano  et  al.  (Juliano  et  al.,  2006). 
They  observed  that  the  elastic  modulus  does  not 
significantly  change  from  PS- 15  to  PS-30: 
(7  /GDe ,«  =  1.17  •  In  this  case,  we  can  assume  that  the 

gel  is  mostly  composed  of  triblock  (80  %  by  mass)  and 
evaluate  x  that  equals  5.3  and  13  for  PS- 15  and  PS-30. 
Assuming  that  changes  in  vA  and  yM  are  negligible  and 

taking  value  of  ^(x,c)  from  our  simulation,  we  can 
estimate  from  Eq.  15  that  GPS_30  /  GPS_l5  =  1.02  • 

These  two  examples  confirm  that  our  approach  can 
predict,  with  reasonable  accuracy,  the  relative  change  of 
elastic  modulus  of  a  micellar  gel  composed  of  triblock 
copolymer  of  high  molecular  weight  with  changing 
polymer  concentration,  relative  length  of  the  blocks, 
solvent  parameters,  and  temperature. 

4.  CONCLUSIONS 


Eq.  14  is  valid  only  for  block  copolymer  of  relatively 
low  molecular  weight  where  entanglement  does  not  likely 
occur.  For  copolymers  with  higher  molecular  weight, 
entanglements  and  attractive  interactions  inside  the 
micelles  would  give  rise  to  a  large  additional  contribution 
to  the  elastic  modulus  and  the  value  of  G  would  be  an 
order  of  magnitude  higher  than  one  predicted  with  Eq. 
(14).  To  address  this  situation  we  propose  to  use  a 
coefficient  depending  on  the  molecular  weight,  Mw  and 
chemical  composition,  q  of  triblock 

copolymer:  p(Mw,q) .  For  this  case,  Eq.  14  becomes 


G  = 


p(Mw,q)RT 

va(x  +  2) 


apb{x,c,%M) 


(15) 


where p(Mw,q)  has  to  be  evaluated  for  each  kind  of 

ABA  block  copolymer  by  comparing  simulation  and 
experimental  results. 

To  verify  this  approach,  we  have  considered  an 
example  of  a  gel  composed  of  SIS  with  higher  molecular 
weight  (Mw=  260  000)  (Laurer  et  al.,  1999).  According  to 
Eq.  15,  with  all  parameters  constant  and  only  changing 
the  volume  fraction  of  SIS  from  0.16  to  0.33  (micellar 
regime)  will  make  the  equilibrium  elastic  modulus  should 
double.  In  Laurer’ s  paper,  G  increases  from  4000  to  9000 
Pa  making  G  2.25  times  higher  for  c  changes  from  0.16  to 
0.33  (density  of  styrene  and  isoprene  blocks  are  taken  to 
be  1.04  and  0.913  g/cm3,  respectively). 


The  equilibrium  and  non-equilibrium  oscillatory- 
shear  DPD  simulation  methods  have  been  used  for 
mesoscale  simulation  of  gel-forming  ABA  triblock 
copolymer  in  a  midblock-selective  solvent.  We  have 
validated  our  approach  using  experimental  data 
demonstrating  that  the  result  of  our  simulation  results  are 
in  good  qualitative  agreement  with  the  experimental  data. 
We  have  developed  a  model  that  estimates  the  absolute 
value  of  the  elastic  response  of  ABA  triblock  copolymer 
of  low  molecular  weight  in  a  B-block  selective  solvent 
with  high  segregation  regime  (%M  >50)  directly  from 
equilibrium  DPD  simulation.  This  approach  eliminates 
the  necessity  of  computationally-costly  non-equilibrium 
shear  simulations. 

We  have  demonstrated  that  our  model  can  predict, 
with  reasonable  accuracy,  the  relative  change  of  the 
elastic  modulus  of  a  gel  consisting  of  triblock  copolymer 
with  changing  polymer  concentration,  relative  length  of 
the  blocks,  solvent  parameters,  and  temperature.  This 
relationship,  derived  from  the  mesoscale  simulation,  will 
help  us  to  reduce  the  time  and  effort  necessary  to  set  up 
experiments  and  eventually  to  control  mechanical 
properties  of  ABA  triblock  copolymer  gels  in  order  to 
mimic  a  wide  range  of  soft  tissues. 

According  to  our  results,  we  can  predict  that  the 
equilibrium  elastic  modulus  of  ABA  triblock  copolymer 
of  low  molecular  weight  would  be  at  its  low  range  for  the 
relative  short  midblock  compared  with  the  endblock 
(more  than  45  %  by  mass  of  styrene  in  dry  copolymer  that 
corresponds  to  x  <  3)  for  any  volume  fraction  of 
copolymer,  c  less  than  0.5.  For  the  lower  styrene  content 
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the  elastic  modulus  does  not  significantly  change  with 
relative  block  size,  and  weakly  decreases  with  decreasing 
styrene  content  after  a  maximum  around  25  %  (x  =  7)  of 
styrene  in  dry  copolymer  for  the  low  copolymer  content 
(c  <  0.2).  For  the  higher  copolymer  concentration  ( c  > 
0.2)  the  elastic  modulus  has  a  pronounced  maximum 
around  25  %  of  styrene  in  dry  copolymer. 

All  of  these  observations  are  valid  for  a  highly 
selective  solvent  with  respect  to  the  midblock.  For  the 
endblocks  with  low  degree  of  polymerization  in  non- 
highly  selective  solvent,  the  equilibrium  elastic  modulus 
has  its  maximum  shifted  to  higher  content  of  polysterene 
and  its  value  drops  significantly  when  the  degree  of 
incompatibility  jM  is  less  than  20.  The  position  of  this 
maximum  should  be  evaluated  for  each  particular  solvent 
from  the  equilibrium  DPD  simulation. 

According  to  our  predictions  we  may  extend  our 
observations  to  the  ABA  copolymers  with  higher 
molecular  weight  including  a  large  additional  contribution 
of  entanglements  and  attractive  interactions  inside  the 
micelles  to  the  equilibrium  elastic  modulus.  However,  to 
predict  the  rheology  of  copolymers  with  higher  molecular 
weight  more  accurately,  we  propose  to  allow 
entanglements  to  occur  in  the  DPD  simulation  (Goujon  et 
al.,  2008)  or  to  exploit  the  “Sliplink”  models  for  block 
copolymers  (Masubuchi  et  al.,  2006). 
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